function [ H_normalized, HPs, beta, gamma, lambda, obj ] = dropoutHMKC(Ks, lambda, nCluster, mid_layer_clusters, mid_layer_width, mid_layer_dropout_ratio)
%*************************************************************************
%
%*************************************************************************
nKernel = size(Ks, 3);
nSmp = size(Ks, 1);
nLayer = length(mid_layer_clusters);
if ~exist('nRepeat', 'var')
    nRepeat = 1;
end
%*************************************************************************
% build layer-by-layer forward connections via layer-wise-dropout
%*************************************************************************
layer_connections = cell(1, nLayer);
for i1 = 1:nLayer
    if i1==1
        n1 = nKernel;
        n2 = mid_layer_width(i1);
    else
        n1 = mid_layer_width(i1-1);
        n2 = mid_layer_width(i1);
    end
    layer_connections{1, i1} = layerDropout(n1, n2, mid_layer_dropout_ratio);
end
%*************************************************************************
% build layer-by-layer backward connections
%*************************************************************************
layer_connections_backward = build_backward_connections(layer_connections, nKernel, mid_layer_width);


%***************************************************************************
%   Repeat 10
%***************************************************************************

best_H = [];
for i_repeat = 1 : nRepeat
    
    HPs = cell(1, nLayer);
    gamma = cell(1, nLayer);
    eta = cell(1, nLayer);
    for i1 = 1:nLayer
        %*************************************************************************
        % initialize all embeddings within the hidden layers following the structure
        %*************************************************************************
        current_layer_input = layer_connections{1, i1};
        %     HPs_current = zeros(nSmp, mid_layer_clusters(i1), mid_layer_width(i1));
        HPs_current = cell(1, mid_layer_width(i1));
        for i2 = 1:mid_layer_width(i1)
            current_input_output = current_layer_input{i2, 1};
            tmp = zeros(nSmp);
            if i1 == 1
                for i3 = 1:length(current_input_output)
                    tmp = tmp + (1/length(current_input_output)) * Ks(:, :, current_input_output(i3));
                end
            else
                HPs_tmp =  HPs{1, i1 - 1};
                for i3 = 1:length(current_input_output)
                    tmp = tmp + (1/length(current_input_output)) * (HPs_tmp{1, current_input_output(i3)} * HPs_tmp{1, current_input_output(i3)}');
                end
            end
            tmp = (tmp + tmp')/2;
            HPs_current{1, i2} = my_kernel_kmeans(tmp, mid_layer_clusters(i1));
        end
        HPs{1, i1} = HPs_current;
        
        %*************************************************************************
        % initialize the connection weights following the structure
        %*************************************************************************
        gamma_current = cell(mid_layer_width(i1), 1);
        for i2 = 1:mid_layer_width(i1)
            gamma_current{i2, 1} = ones(1, length(current_layer_input{i2, 1})) / sqrt(length(current_layer_input{i2, 1}));
        end
        gamma{1, i1}= gamma_current;
        
        %*************************************************************************
        % initialize the connection weights following the structure
        %*************************************************************************
        eta{1, i1}= ones(1, mid_layer_width(i1))/sqrt(mid_layer_width(i1));
    end
    %*************************************************************************
    % initialize the final connection weight
    %*************************************************************************
    beta = ones(1, mid_layer_width(end))/sqrt(mid_layer_width(end));
    
    
    iter = 0;
    not_converges = 1;
    objHistory = [];
    while not_converges
        %%
        %*************************************************************************
        % update the final consensus result H
        %*************************************************************************
        H = update_H(beta, HPs{nLayer}, nCluster);
        
        %%
        %*************************************************************************
        % update the hidden embeddings
        %*************************************************************************
        %     HP = update_HP_nor(H, HP, Ks, beta, gamma, lambda, mid_layer_clusters);
        %     HP = update_HP_rev(H, HP, Ks, beta, gamma, lambda, mid_layer_clusters);
        
        if rem(iter, 2) == 0
            HPs = dropout_update_HP_forward(H, HPs, Ks, beta, gamma, eta,  mid_layer_clusters, mid_layer_width, layer_connections, layer_connections_backward);
        else
            HPs = dropout_update_HP_backward(H, HPs, Ks, beta, gamma, eta,  mid_layer_clusters, mid_layer_width, layer_connections, layer_connections_backward);
        end
        
        %%
        %*************************************************************************
        % update the final layer weights
        %*************************************************************************
        beta = update_beta(H, HPs{nLayer});
        
        %%
        %*************************************************************************
        % update the hidden layer weights gamma
        %*************************************************************************
        gamma = update_gamma(HPs, Ks, layer_connections, mid_layer_width);
        
        %%
        %*************************************************************************
        % update the hidden layer weights eta
        %*************************************************************************
        eta = update_eta(HPs, Ks, gamma, layer_connections, mid_layer_width);
        
        %%
        iter = iter+1;
        obj = calculate_obj(H, HPs, Ks, beta, gamma, eta, mid_layer_width, layer_connections);
        objHistory = [objHistory; obj]; %#ok
        
        if (iter >= 10) && (abs( (objHistory(iter - 1) - objHistory(iter)) /objHistory(iter) ) < 1e-5 || iter > 100)
            not_converges = 0;
        end
    end
    
    if isempty(best_H)
        best_H = H;
        best_HPs = HPs;
        best_gamma = gamma;
        best_eta = eta;
        best_beta = beta;
        best_lambda = lambda;
        best_obj = objHistory;
        best_objend = calculate_obj(H, HPs, Ks, beta, gamma, eta, mid_layer_width, layer_connections);
    else
        objend = calculate_obj(H, HPs, Ks, beta, gamma, eta, mid_layer_width, layer_connections);
        if best_objend > objend
            best_H = H;
            best_HPs = HPs;
            best_gamma = gamma;
            best_eta = eta;
            best_beta = beta;
            best_lambda = lambda;
            best_obj = objHistory;
            best_objend = calculate_obj(H, HPs, Ks, beta, gamma, eta, mid_layer_width, layer_connections);
        end
    end
    
end
H = best_H;
HPs = best_HPs;
gamma = best_gamma;
eta = best_eta;
beta = best_beta;
lambda = best_lambda;
obj = best_obj;
H_normalized = H;
end

